A pulsating regime of magnetic deflagration 
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The stability of a magnetic deflagration front in a collection of molecular magnets, such as Mni2-acetate, is 
considered. It is demonstrated that stationary deflagration is unstable with respect to one-dimensional perturba- 
tions if the energy barrier of the magnets is sufficiently high in comparison with the release of Zeeman energy 
at the front; their ratio may be interpreted as an analogue to the Zeldovich number, as found in problems of 
combustion. When the Zeldovich number exceeds a certain critical value, a stationary deflagration front be- 
comes unstable and propagates in a pulsating regime. Analytical estimates for the critical Zeldovich number are 
obtained. The linear stage of the instability is investigated numerically by solving the eigenvalue problem. The 
nonlinear stage is studied using direct numerical simulations. The parameter domain required for experimental 
observations of the pulsating regime is discussed. 



I. INTRODUCTION 



Molecular magnets have been a subject of intense exper- 
imental and theoretical studies for almost two decades, see 
Refs. OH for a review. The field belongs to mesoscopic 
physics as the typical length scales are between micro- and 
macroscopic. Such an intermediate parameter regime pro- 
vides unique conditions for observing quantum and classi- 
cal phenomena acting together. Besides, a large magnetic 
spin of a single molecule makes molecular magnets natural 
candidates for novel magnetic storage media and quantum 
computing^. 

Two prominent and well investigated materials in this re- 
search field are Mn^-acetate and Fes, both which demonstrate 
super-paramagnetic properties 1 . They both have a large spin 
at the ground state (S = 10) and strong magnetic anisotropy 
with 10 pairs of degenerated levels corresponding to positive 
and negative projections S z on a chosen axis, and an addi- 
tional level with S z = 0^"—. When an external magnetic field 
is applied to a sample of molecular magnets along the crystal 
axis, then spin orientation in the direction of the field becomes 
preferable. At low temperature all molecules populate only 
one level (e.g. S z = 10) and the magnetization reaches the sat- 
uration value. If the direction of the external magnetic field 
changes to the opposite one, then the former ground state of 
the molecules becomes metastable with an increased potential 
energy (the Zeeman energy). An energy barrier hinders direct 
transition of the molecules from the metastable state to the 
new ground state, which, therefore, can occur as a thermal re- 
laxation or as avalanches. Thermal relaxation goes slowly and 
uniformly in space for the whole sample. At temperatures of 
a few Kelvin this process may be neglected, since its charac- 
teristic time is rather large, e.g., about two months for Mni2 at 
2 K 5 . However, many experimental works on molecular mag- 
nets demonstrated quite fast transition in a form of avalanche 
with the characteristic time about few milliseconds^^. Re- 
cent detailed studies of the avalanche revealed that the spin re- 
versal does not happen simultaneously within the whole sam- 
ple in that case, but it occurs in a narrow zone propagating as 
a front at a velocity of several meters per seconcU-L - — . The 
avalanche is accompanied by significant heat release similar 
to a deflagration wave in combustion and for this reason it 
was named "magnetic deflagration". 



In slow combustion, the deflagration front corresponds to a 
thin zone of chemical reactions separating the cold fuel mix- 
ture and the hot burned products 19 . The released energy is 
transmitted to the cold fuel mixture due to thermal conduc- 
tion; the temperature of the fuel mixture increases, which 
stimulates chemical reactions. The slow combustion front 
propagates in gaseous mixtures with a substantially subsonic 
velocity within the range from several centimeters to sev- 
eral meters per second 19 . Analogously, in magnetic deflagra- 
tion, the energy of chemical reactions is replaced by the Zee- 
man energy of the molecular magnets in an external magnetic 
field. The released energy is distributed to the neighboring 
layers by thermal conduction, which increases the tempera- 
ture of the originally cold medium and leads to much higher 
probability of spin transition. In recent studies, Garanin and 
Chudnovsky 12 developed a theoretical description of a pla- 
nar stationary magnetic deflagration with maximal possible 
release of the Zeeman energy (they called such a regime "full 
burning"). Since then, a number of papers have been devoted 
to ignition techniques and speed measurements of the mag- 
netic deflagration in Mni^"^. Magnetic avalanches have 
been also observed in other materials, like the intermetal- 
lic compound GdgGe/p22. However, the magnetic "burning" 
does not necessarily have to be complete. As we show in the 
present paper, the Zeeman energy release depends on the ex- 
ternal magnetic field and on the initial concentration of active 
molecules. A sufficiently low energy release may lead to a 
new pulsating regime of magnetic deflagration. 

The aim of the present work is to study the stability of a 
planar stationary front of magnetic deflagration. Here, we 
demonstrate that stationary deflagration is unstable with re- 
spect to one-dimensional perturbations if the Zeeman energy 
release at the front is sufficiently low in comparison with the 
energy barrier of the spin transition. The condition of the in- 
stability may be formulated using the Zeldovich number, sim- 
ilar to the case of combustion of solid propellant a l9 i 21 ~ 2sL. In 
the problem of magnetic deflagration, the Zeldovich number 
is essentially given by the ratio of the energy barrier of the 
molecular magnets and the temperature at the deflagration as 
determined by the Zeeman energy release. When the Zel- 
dovich number exceeds a certain critical value, a stationary 
deflagration front becomes unstable and propagates in a pul- 
sating regime. We obtain analytical estimates for the critical 
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FIG. 1 : The energy levels for the magnetic molecule system of Mn]2 
in the external magnetic field H z = 0.2T 



Zeldovich number. We investigate the linear stage of the in- 
stability numerically by solving the eigenvalue problem. We 
study the nonlinear stage using direct numerical simulations. 
We also discuss the experimental parameters required for ob- 
servations of the pulsating regime. 



II. A PLANAR STATIONARY FRONT WITH 
INCOMPLETE ZEEMAN ENERGY RELEASE 

We consider a system of molecular magnets placed in an 
external magnetic field taking Mni2-acetate as a particular ex- 
ample. The simplified Hamiltonian of the system has been 
suggested to take the formic 



3tf = -DS%- g\i B H z S z 



(1) 



where S denotes the spin, D sa Q.65K is the magnetic 
anisotropy constant, g w 1 .94 is the gyromagnetic factor, /ig is 
the Bohr magneton, and H z is the external magnetic field. The 
energy levels for a molecule magnet Mni2 in the external field 
of H z = Q.2T are depicted in Fig. 1. In Fig. 1 we indicate the 
Zeeman energy Q and the energy barrier E u of the transition 
from the metastable state to the ground state, which plays the 
role of the activation energy for magnetic deflagration (both 
values are presented in temperature units). Using the Hamil- 
tonian Eq. (Q]i we find how these two energies depend on the 
magnetic field, according to 



and 



E a = DS 2 l0 - gH B H z S w + j^V-ltf 



Q = 2gll B H z S z . 



(2) 



(3) 



The last term in Eq. d2j is quite small as compared to the 
first two terms, so the activation energy decreases with in- 
crease of the magnetic field. However, the Zeeman energy 
increases linearly with the magnetic field. If the magnetic 



field is high enough, the energy difference between the ground 
state (S z = — 10) and the metastable state (S z = 10) is rather 
large, so that all the molecules tend to occupy the lowest en- 
ergy level, i.e. S z = — 10 in Fig. 1. When the direction of the 
external magnetic field switches to the opposite one, then the 
ground state and the metastable state exchange places and all 
the molecules tend to change their spin projection to the oppo- 
site one. Garanin and Chudnovsky identified such a transition 
as "full burning" of molecular magnets 12 . They also indicated 
that the regime of full burning is possible if the magnetic field 
is stronger than a certain critical value. Still, the external mag- 
netic field is a controlled parameter of the experiments, which 
allows the possibility of incomplete burning in magnetic de- 
flagration. We also point out that the analytical method used 
to calculate the magnetic deflagration speed in 1.1211 is rigorous 
only in the limit of an activation energy large compared to the 
energy release, a condition which is satisfied to a much better 
degree in the regime of incomplete burning. In the present pa- 
per we are interested in the regime of incomplete burning in 
magnetic deflagration achieved for a sufficiently low external 
magnetic field. In that case the Zeeman energy is rather low, 
so that we obtain the fraction of molecules on the first energy 
level according to 



1 



ex P (e/r) + r 



(4) 



where total number of molecules in all energy levels corre- 
sponds to unity. Still, in the typical experimental conditions, 
population of levels above the first one is negligible. At the 
same time, the fraction at the first level, Eq. ©, may be sig- 
nificant in a low magnetic field and cannot be neglected. As 
we will see, the stability of magnetic deflagration is sensitive 
to the activation energy scaled by the energy release in the pro- 
cess. This ratio increases also when the active Mni2 molecules 
are "diluted" with some neutral media as it was performed ex- 
perimentally in, e.g., ll26l - r28ll . When mixed carefully with 
other compounds, the Mnj2 molecules retain their magnetic 
properties. Since the heat release in magnetic deflagration 
happens only due to active magnet molecules, then the neu- 
tral compound reduces total energy release and the final tem- 
perature of the sample, Tf, thus increasing the scaled activa- 
tion energy E a /Tf. As a result, the scaled activation energy in 
magnetic deflagration becomes a free parameter, which may 
be controlled in the experiments by the external magnetic field 
and by the level of dilution. 

The governing equations for magnetic deflagration are 12 



dE 
dt 



V-( K VE)-fQ 



d7' 



dn 



1 



exp 



Ea 
T 



ex P (e/r) + iJ 



(5) 



(6) 



where E is the thermal phonon energy, K is thermal diffu- 
sion of energy, Q is the Zeeman energy release in tempera- 
ture units, / < 1 indicates possible reduction of the energy 
release because of dilution of active molecules (/ = 1 corre- 
sponds to zero dilution), n is fraction of magnetic molecules 



3 



in the metastable state, E a is the energy barrier of tunneling 
measured in temperature units and %r is a constant of time 
dimension. In magnetic deflagration, thermal diffusion and 
heat capacity depend strongly on temperature. Following 12 , 
we take the heat capacity in the classical form corresponding 
to phonons^ 



C = Ak, 



&D 



(7) 



where ©£> is the Debye temperature, with ©o = for Mni2, 
kg is the Boltzmann constant, A = 127T 4 / 5 corresponds to the 
simple crystal model, a is the problem dimension (we take 
a = 3 in most of the calculations, which corresponds to the 3D 
geometry). Dependence of thermal diffusion on temperature 
may be taken in the form K °= T~P, where parameter j3 was 
considered in Refs. Il2ll30ll within the range between and 
13/3. In the present paper we take j3 within the same range, 
and show that it has minor influence on the deflagration sta- 
bility. Using the definition of heat capacity, C = dE/dT, we 
find the phonon energy as a function of temperature 
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FIG. 2: The final temperature and the scaled activation energy in 
magnetic deflagration versus the magnetic field for the initial tem- 
perature Tq = 0.2K and two dilution factors / = 1; 1/3. 
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An important parameter of deflagration dynamics is deter- 
mined by the ratio of the energy barrier (in temperature units) 
and the temperature in the hot region, E a /Tf. A combustion 
counterpart of this value is related to the activation energy of 
chemical reaction, which is typically rather large. In the case 
of complete burning in magnetic deflagration of Mni2 this pa- 
rameter was evaluated in IU2I1 as E a /Tt ~ 6, which is not a 
large value. At the same time, this parameter may increase 
almost without limits at low magnetic fields and for consid- 
erable dilution of active molecular magnets. When the ratio 
E a /Tf is very large, the transition from the metastable to sta- 
ble states goes in a thin region where the Arrhenius function in 
Eq. © is different from zero. In this limit we may obtain an 
analytical solution to Eqs. (0, (O by the classical method of 
Zeldovich-Frank-Kamenetskyi2, which was applied to mag- 
netic deflagration in 111 211 . In a general case of finite E a /Tr, the 
system Eqs. ©, © may be solved numerically as an eigen- 
value problem. 

First, we calculate the final temperature in the hot region 
behind the magnetic deflagration front. Using the energy con- 
servation and Eqs. @, (0, ([8]l we obtain 



Ar a+1 



/g(a + l)0g 



1 



1 



AT, 



a+l 



fQ(a 



exp(Q/T ) + l 
1 

Tj0g + e X p(e/7>) + l' 



(9) 



Equation (0 determines the final temperature Tf as a function 
of the magnetic field, the dilution factor / and other parame- 
ters of the process. The first terms on both sides stand for the 
thermal energy, though initial thermal energy is usually rather 
small. The other terms indicate the fractions of molecules oc- 
cupying the first energy level. In the case of full burning con- 
sidered in 1U2I1 this equation becomes much simpler and may 



be solved analytically. In the case of incomplete burning, a 
numerical solution to Eq. (0 is required. In Fig. 2 we show 
the numerical solution to Eq. (0, i.e. the final temperature 
and the scaled activation energy versus the magnetic field for 
two dilution factors / = 1; 1/3. In Fig. 2 we take the initial 
temperature Tq = Q.2K, though it has minor influence on the 
result. The dilution factor / = 1 stands for pure Mni2 me- 
dia, while / = 1/3 means that the average energy release is 
3 times lower in comparison with the pure substance. The fi- 
nal temperature increases with the magnetic field, while the 
scaled activation energy decreases. The change in the scaled 
activation energy due to the dilution factor may be also quite 
strong. 

We consider a stationary solution to Eqs. (|5}, © in the 
form of a planar front propagating with velocity Uf. To be 
particular, we assume that the front moves along the x-axis in 
the negative direction as shown in Fig. 3. Taking the reference 
frame of the front, we find 



d , d 



dE\ 



TJ — - ~— ( - — 

f dx~~T R eXP \ T 



exp(Q/T) 



1 



(10) 



(11) 



Integrating Eq. dTOb from the initial to final states (labels 
and f, respectively) and neglected the initial thermal energy, 
we obtain the final energy as 



Ef=fQ(n -n f ) 



(12) 



Integral to Eq. ( fTOb specifies the internal structure of the tran- 
sition zone 



. s dE 

UffQ(n Q -n f )=K f — . 



(13) 
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In the limit of a large activation energy, E a /Tr^> 1, the transi- 
tion region may be presented as a surface of weak discontinu- 
ity, where energy and temperature are continuous and tend to 
their maximal values E — > Ef, T — » Tf(Ef), while their deriva- 
tives experience jump. In this limit one obtains the magnetic 
deflagration velocity^ 



I K f ( E c 



2T, 



where 



E a fQ(n f -no) 



C f T, 



(a + 1) 7) 



(14) 



(15) 



plays the role of the Zeldovich number; the last relation in 
Eq. CUl follows from Eqs. ©, ® and (H2J. Equation (Q3]i 
for the deflagration velocity coincides with Eq. (109) in Ref. 
Illl . With the accuracy of a factor of 1/ (a + 1), the Zel- 
dovich number shows the ratio of the activation energy and 
the final temperature in the deflagration front. Strictly speak- 
ing, the approach of a thin transition zone holds as long as 
Zeldovich number is large Z> 1, which is not applicable to 
the regime of full burning. We can find the dependence of the 
Zeldovich number upon the main experimental parameters of 
the problem taking into account Eqs. to,©,©,® and ( TT2| > 



2(a + 1 ) a + 2 %gli B S z H z f (n Q - n f ) 



DS 2 -g^ B H z S w + ^ 2 H 2 ). (16) 



Equation (IT6b shows how the Zeldovich number depends on 
the magnetic field and the dilution factor. In particular, high 
values of the Zeldovich number correspond to the low field 
strength. For a simplified quantitative estimate one can ne- 
glect the magnetic terms in the second couple of parentheses, 
which provide a contribution less than few percents for fields 
about 0. 1 T and smaller, and find an evaluation for the Zel- 
dovich number 



vary because of the front perturbations, and the instant front 
velocity may be written as 



U, = t/yexp 



2T, 



f 



Ea_ 

2T, 



(18) 



where label / refers to the stationary case, whilst t indicates 
the time-dependent temperature in the infinitely thin transition 
zone. Position of the transition zone, x — 0(f), in the station- 
ary reference frame is determined by the equation 



d(j) 
It 



~57=- U f 



2T, 



f 



IT, 



1 



(19) 



where the first minus sign indicates propagation of the front in 
the negative direction. Boundary conditions at the transition 
region may be found similar tc 19 - 22 . Integrating Eq. (O twice 
over the transitional zone, we obtain continuous energy and 
temperature 



Ei,-L. — &)_, 



(20) 



where the labels (j)± correspond to hot and cold sides of the 
transition zone. The first derivative of energy (i.e. energy 
flux) experiences jump at the interface as 



U f Qn exp 



fdE 



2T f 



Ea_ 

2T, 
fdE 



(21) 



We investigate the linear stability of the stationary deflagration 
and consider small perturbations of all variables in the expo- 
nential form E oc exp(<7f ), where the growth rate a may have 
both real and imaginary parts. We solve the stability prob- 
lem analytically in the limit of a thin transition zone similar 
to 1 19]. Outside the transition zone the perturbed equations 
©-dUl are 



dE d 2 , ~\ 



(22) 



Z^l.3[HJ (»o-«/)] 
where H z is the magnetic field in tesla. 



-1/4 



(17) 



III. ANALYTICAL ESTIMATES FOR THE PULSATION 
INSTABILITY 

In this section we obtain an analytical scaling for the ID 
instability of the magnetic deflagration front in the model 
of a thin transition zone, i.e. at a large Zeldovich number 
Z> 1. In that case, according to Eq. (14\ . the deflagra- 
tion velocity is extremely sensitive to temperature variations 
in the transition zone. In comparison with the Arrhenius func- 
tion exp (—E a j2Tf), all other parameters in Eq. (Q3) may be 
treated as constant. Temperature in the transition zone may 



dn 

an + Uf— = 0. 
ox 



(23) 



Behind the transition zone we have n ~ 0. Similar to the 
respective combustion problem^, we consider first a hypo- 
thetic case of constant coefficient of thermal conduction K = 
const = Kf. Solving Eq. (l22l outside the transition zone with 
E oc exp(jj,x) we find two modes 



2 U f G n 

M -— M =0, 

*f K f 



JU0,1 



El 

2k f 



± 
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II 2 
4k 2 



a 



(24) 



(25) 



5 



with labels and 1 indicating media ahead and behind the 
transition zone. In the case of magnetic deflagration, thermal 
conduction decreases strongly with temperature, K oc T~P. 
Still, the mode behind the transition zone with ji\ < is an ex- 
act solution to the linearized Eq. ([T9T > even in that case, since 
temperature behind the transition zone is uniform in the sta- 
tionary deflagration. On the contrary, the mode ahead of the 
transition zone with /io > is only an approximation. Simi- 
lar to the combustion theory 19 , we may define the parameter 
Lf= Kf/Uf as the deflagration front thickness related to ther- 
mal conduction in the hot region. At the same time, the char- 
acteristic length scale of temperature profile in the stationary 
deflagration increases in the cold layers as k(T)/U/, which 
allows to consider them as quasi-uniform with respect to the 
perturbations. Indeed, the mode ahead of the transition zone 
describes decay of perturbations at the length scale Hq with 



1 

7L, 



1 + Jl+AoL f /U f 



7^- (26) 



Therefore, investigating the ID instability analytically, it is 
justified to consider only the heating layer of the size about 
Lf close to the transition zone and to treat the coefficient of 
thermal conduction as approximately constant. The numerical 
solution to the problem obtained below supports the analytical 
approximation. 

We consider perturbations of the boundary conditions (|T9l 
- (ETb at the transition surface 



as a function of the Zeldovich number. According to Eq. (|3T1 i. 
the instability develops for Z > 4 + 2\/5 ~ 8.5. Close to this 
critical value the real part of the growth rate Red goes to zero, 
while the imaginary part remains finite, co = Ima ^ 0, which 
indicates the pulsation regime of the instability. The obtained 
critical value of the Zeldovich number corresponds to rather 
high scaled activation energy, E a /Tf w 34, see Eq. (1151 , in 
comparison with the values typical for the regime of full burn- 
ing Ref. lfl2tl . Still, this value is attainable experimentally 
for smaller magnetic fields and some dilution of the active 
molecules. At even larger value of the Zeldovich number, 
Z = 11.7, corresponding to E a /Tf = 46.8, Eq. (fJTJ demon- 
strates bifurcation, so that the instability growth rate becomes 
purely real with zero imaginary part for Z > 1 1 .7. Below we 
will compare the analytical scaling Eq. (|3T| i to the numerical 
solution to the problem taking into account finite width of the 
transition zone. 



IV. NUMERICAL SOLUTION TO THE STABILITY 
PROBLEM TAKING INTO ACCOUNT FINITE WIDTH OF 
THE TRANSITION ZONE 

In this section we solve the stability problem numerically 
taking into account finite width of the transition zone. We 
introduce dimensionless variables and parameters 

6 = T/Tf, a = n/n , %=x/L f , 



-u, 



2C f Tj 



f dE 



dE 

x 



E a ( dE _ 



J dlE 



The last equation may be also modified as 

E a - ~ ~ Q 

- QnQ——rEp + = HiL f E^ + - ^L f E^ - <j)n Q — 

2tfl f Lf 



(27) 



E^,+ = E^ + (j) 1'— ) =E^ + <j)noQ/L f . (28) 



(29) 



taking into account the perturbation modes. Thus, we have 
to solve the algebraic system of three equations (l27l . (l28t 
and ( f30b , taking into account the modes d25l ). After heavy 
but straightforward algebra we end up with a simple quadratic 
equation 



—J - — (Z 2 /4-2Z-l)+2Z = 0, (31) 



which is similar to the respective result in the combustion 
theory 19 . Equation OTT ) describes the instability growth rate 



© = Ea/Tf, A = Q/T f , K=K f 9- 
and rewrite the evolution equations (E)-© to 

1 



(32) 



R ad9 a dd _ d 
6 d^ + d d$~M 



+ 



./Acxp j -® 



exp(A/0) 



1 



da dc 



_ + _ = _ Aexp ^__ 







1 



exp(A/0) 



(33) 



(34) 



where A — Lf/ (trUA is an eigenvalue of the stationary prob- 
lem and the designation J corresponds to the ratio of Zeeman 
and thermal energies, 



fQ% 

Ak B T f a+l 



(35) 



It can be shown that the parameter J is almost a constant, 
/~l/(a + l). We also introduce the thermal flux, y = 
9 a ~^d9/d^, so that the stationary deflagration is described 
by the system of equations 



3^- = p \j/- JAa exp - — 



(36) 
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FIG. 3: Stationary profiles of concentration, temperature and energy 
release. The plot parameters are E a /Tf = 30, /3 = 3, Sq = 0.2 



Typical profiles of scaled concentration, temperature and en- 
ergy release 



/ © 

W = Aexp I — — 



fl - ( ex P( q I + 1 



(37) 



are depicted in Fig. 3 for the magnetic deflagration propagat- 
ing to the left with j3 =3. In Fig. 3, the scaled activation en- 
ergy is taken rather high, © = 30, with the initial temperature 
00 = 0.2. The front velocity is determined by the eigenvalue 
of the problem A in Eq. ( 1361 ). which was computed numeri- 
cally using the shooting method similar to2i. We can see in 
Fig. 3 that final number of molecules in the metastable state 
behind the deflagration front is different from zero «/ = 0.47, 
which indicates that burning is incomplete. Initial number of 
molecules in the metastable state is different from unity too, 
«o = 0.77, due to the non-zero temperature ahead of the front 
before switching of the magnetic field, see Eq. (0). Figure 
3 shows also that the total thickness of the deflagration front 
is much larger than Lf. This difference in the characteristic 
length scales should be attributed, first of all, to the temper- 
ature dependence of thermal conduction. Still, even in the 
case of constant thermal conduction typically used in combus- 
tion problems, the effective flame thickness is almost order-of- 
magnitude larger than the conventional definition for Lf, e.g. 
sseM. The asymptotic temperature behavior ahead of the front 
is described by a single differential equation, derived from Eq. 
(l33l , assuming that exp(— 0/0) << 1 in the cold matter 



de_ 



a + l 



< +1 



9 a+l ) =0. 



(38) 



Equation ( 1381 1 determines variations of the length scale in the 
temperature profile. 

Now we consider stability of the magnetic deflagration. We 
investigate dynamics of small perturbations of the stationary 
solution, so that all variables may be presented as <p = <f>o + 



<pexp (yr). Then the perturbed system (l36l l is 
<9y/ 

de_ 

W 



= 9 p y/ + (9 a a+p9^ 1 \ir-JG)9-JAsxp 
= e^^-r- (j3 - a) e^^ye 



where the designation 



(39) 



G = Aexp 



9)9 
& 



2 [®a- 



exp(A/0) + l / exp(A/0) + l 



has been introduced for brevity. Parameter y stands for the 
scaled instability growth rate, y = aLf/Uf. In order to spec- 
ify boundary conditions for the numerical solution, we find 
the decaying perturbation modes in the uniform regions, <p °c 
(f)exp(p^). The system (1391 involves three modes in the uni- 
form regions: one in the hot matter with fj, < 0, % — » +°°, 
and two in the cold matter with jj. > 0, £, — s- — °°. We inte- 
grate the system d39l numerically three times corresponding 
to the three perturbation modes: two times from the left (cold) 
side and one from the right (hot) side. At some point, e.g., at 
the maximum of the energy release, the solutions form a ma- 
trix, with the determinant depending on the instability growth 
rate, y. The condition of zero matrix determinant specifies the 
growth rate y as an eigenvalue of the system d39l . This method 
has been applied successfully in studies of different hydrody- 
namic instabilities, e.g. the Rayleigh-Taylor instability and 
the Darrieus-Landau instability of combustion and laser abla- 
tion, as well as for other plasma instabilities 31 -^"—. 



V. RESULTS AND DISCUSSIONS 

Numerical solution to the stability problem is shown in Fig. 
4 for different values of the scaled activation energy propor- 
tional to Zeldovich number as E a /Tf = 4Z for a 3D problem; 
other deflagration parameters are the same as used in Fig. 3. 
As we can see in Fig. 4, planar stationary magnetic deflagra- 
tion is unstable at sufficiently large values of the scaled activa- 
tion energy; the stability limit was calculated as E a /Tr — 28.2. 
The instability domain consists of two parts separated by the 
bifurcation point at E a /Tf — 33.2. Most of the domain (to the 
right of the bifurcation point) corresponds to purely real in- 
stability growth rate a > with two branches describing the 
fast and slow perturbation modes (shown by the solid lines). 
The instability growth rate of the fast mode increases with 
the scaled activation energy without any limit. The analyt- 
ical model, Eq. QTb , predicts the asymptotic increase of the 
growth rate for the fast mode as ff« (E a /Tf) 2 for E a /Tt — 
The growth rate of the slow mode decreases with the scaled 
activation energy. For this reason it is expected that the slow 
mode plays a noticeable role only close to the bifurcation 
point, where the growth rates of the fast and slow modes are 
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FIG. 4: The instability growth rate cr and the perturbation frequency 
ft) versus the scaled activation energy calculated for the same param- 
eters as in Fig. 3. Solid lines show the domain of zero frequency; the 
dashed and dotted lines correspond to the regimes when perturba- 
tions have both real (growth rate, a) and imaginary part (frequency, 
ft)) in the problem eigenvalue. The dash-dotted lines present the an- 
alytical model Eq. l !28t . The markers show results of the direct nu- 
merical simulations: empty circles stand for stable region and filled 
triangles represent the unstable pulsating regime of the deflagration. 



comparable. In a small part of the instability domain, for inter- 
mediate values of the scaled activation energy in between the 
stability limit and the bifurcation point, 28.2 < E a /Tf < 33.2, 
the instability growth rate is complex with the real part, a > 
(dashed line), and imaginary part, co (dotted line). Although 
this part of the instability domain is rather small, it indicates 
the physical outcome of the instability at the nonlinear stage. 
Because of the non-zero frequency, (0 ^ 0, it is natural to ex- 
pect that the instability leads to a pulsating regime of mag- 
netic deflagration at the nonlinear stage. The analytical solu- 
tion, Eq. OT), obtained within the model of a thin transition 
zone (dashed-dotted lines) shows qualitatively the same stabil- 
ity properties of magnetic deflagration as the numerical solu- 
tion. Still, we observe minor quantitative difference between 
the analytical model and the numerical solution. According 
to the analytical model, the stability limit and the bifurcation 
point are expected at E a /Tf = 34.0 and E a /Tf = 46.8, which 
differ by approximately 20-30% from the respective numeri- 
cal values. The limited accuracy of the analytical model is due 
to the simplifying approximations of 1 ) a constant coefficient 
of energy diffusion and 2) an infinitely thin zone of energy 
release. In particular, the shortcomings of the discontinuity 
model have been discussed in |24||25I1 in the context of solid 
propellant combustion. 

The numerical solution to the stability problem indicates 
the experimental parameters required to observe the unstable 
non-stationary regime of magnetic deflagration. We plot the 
stability diagram in Fig. 5 in / — H z coordinates using Eq. 
( fTTI i; the shading represents absolute value of the instability 
growth rate. The dashed white line in Fig. 5 corresponds to 
the critical value of the Zeldovich number (the stability limit) 




0.0001 0.001 0.01 0.1 



H(T) 

FIG. 5: The stability limit (dashed line) of the magnetic deflagration 
in coordinates of the dilution factor / and the magnetic field H. The 
shading shows the instability growth rate. Other solid lines corre- 
spond to the respective constant values of the Zeldovich number. 



obtained numerically. Magnetic deflagration propagates sta- 
tionary in the parameter domain to the right of the stability 
limit, in the region of a high magnetic field. To the left of the 
critical curve, for a low magnetic field, the stationary magnetic 
deflagration is unstable, and we expect a pulsating regime of 
the deflagration front. Particularly, in the case of Mni2 with 
the dilution level / = 1/3 one should expect the instability for 
the magnetic fields below 10~ 2 T, which is possible to achieve 
experimentally. 

In order to understand front dynamics at the nonlinear stage 
of the instability, we perform direct numerical simulation of 
Eqs. d33l-(l34ll. In our simulations we use the method of finite 
difference for the spatial derivatives and the common fourth 
order Runge-Kutta method for the time step integration. Since 
the front structure is ID, then we were able to obtain fine front 
resolution together within acceptable time of computational 
runs. The results of our simulations, i.e. evolution of the 
magnetic deflagration speed, are shown in Fig. 6 for different 
values of the scaled activation energy. Figures 6 a, b demon- 
strate front behavior close to the stability limit in the stable 
(a) and unstable (b) parameter domain. In the first case with 
E a /Tf = 28.55, the velocity perturbation oscillates and decays 
in time, which implies a negative instability growth rate. On 
the contrary, in Fig. 6b with E a /Tf = 28.8, the amplitude 
of velocity oscillations grows in time, which corresponds to 
the instability onset. The markers in Fig. 4 show the oscilla- 
tion frequency of the magnetic deflagration obtained in direct 
numerical simulations with empty circles and filled triangles 
standing for the stable and unstable regimes, respectively. As 
we can see, the direct numerical simulations are in a good 
agreement with the solution to the eigenvalue problem, which 
concerns both the stability limits and the oscillation frequency. 
Still, the instability is rather weak in Fig. 6 b with the os- 
cillations described well by the sine function. The nonlinear 
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effects become noticeable for higher values of the scaled ac- 
tivation energy with sharp peaks and smooth troughs in the 
oscillations as presented in Fig. 6c for E a /Tf = 30.2. The 
simulations explain the physical mechanism of the obtained 
instability. Magnetic deflagration propagates due to two ef- 
fects: release of the Zeeman energy in a thin transition zone 
and transfer of the energy to the cold layers (preheating) due 
to thermal conduction. In a stationary regime of magnetic de- 
flagration, these two processes work at the same rate and bal- 
ance each other. However, at high values of the activation 
energy, the rate of energy release is too sensitive to temper- 
ature in the transition zone. Small temperature perturbations 
may increase the "burning" rate considerably, which makes 
the transition zone sweep fast over the preheated matter until 
it comes to cold regions and stops waiting for a new portion 
of cold matter to be preheated. As soon as it happens, next 
pulsation of the front takes place. 

It is interesting that the three figures 6 a, b and c demon- 
strate three qualitatively different regimes of magnetic defla- 
gration, though the activation energy changes only within 5% 
from plot (a) to plot (c). In Fig 6 d we take the activation 
energy E a /Tf = 33.8 close to the bifurcation point and see a 
complicated front behavior at the onset of the velocity pulsa- 
tions. Presumably, the complicated behavior happens because 
of two unstable modes with close instability growth rates at 
the vicinity of the bifurcation point. Still, after an initial tran- 
sition period, the front pulsations resemble those of Fig. 6 
c. Finally, taking the activation energy noticeably larger than 
the bifurcation point, E a /Tf = 40, we observe even more pro- 
nounced nonlinear features in the front oscillations with even 
sharper peaks of large amplitude, shown in Fig. 6e. The ob- 
tained results are qualitatively similar to the pulsation insta- 
bility of solid propellant combustion studied in, e.g., ll2lll23ll . 

As we have shown, the scaled activation energy (or the Zel- 
dovich number) is the main parameter, which controls ID sta- 
bility properties of magnetic deflagration. Still, the problem 
involves other parameters as well, namely, power exponents 
a and j3 in the heat capacity and thermal conduction, respec- 
tively, and the scaled initial temperature of the system 6q. Fig- 
ure 7 shows the instability growth rate for different types of 
thermal conduction (different values of the power exponent 
j8, 13/3>/3>0). As we can see, a particular type of the 
thermal conduction influences the temperature profile in the 
stationary magnetic deflagration noticeably in agreement with 
Eq. (1341 , see the insert of Fig. 7. At the same time, the sta- 
tionary profiles of concentration and energy release remain the 
same for different values of the j3 -factor. The power exponent 
j3 determines the total length scale of the deflagration front in 
the cold region. However, even considerable variations of j3 
modify the instability growth rate only slightly: the stability 
limit changes only by 5% from E a /Tf = 27.9 to E a /Tf = 29 A; 
the bifurcation point changes within 15% from E a /Tf = 32 to 
Ea/Tf = 37. This result demonstrates that the perturbation 
modes are strongly localized around the zone of energy re- 
lease within the length scale about Lf as we explained in Sec. 
III. For this reason, considerable modifications of the temper- 
ature profile in the far zone on the length scales much larger 
than Lf due to variations of j3 produce a minor effect on the 
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FIG. 6: Deflagration speed versus time for different values of the 
scaled activation energy E a /T f = 28.55; 28.8; 30.2; 33.8; 39.6 for 
plots (a) - (e), respectively. 
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FIG. 7: The instability growth rate o" versus the scaled activation 
energy calculated for different values of the power exponent |3. Other 
parameters are a = 3, 8q = 0.2. The inset presents the respective 
temperature profiles for E a /Tj = 30 and /3 = 0; 3. 



stability properties. Influence of the initial temperature on the 
instability development is even weaker. In particular, chang- 
ing the scaled initial temperature from 6q = 0.2 to 9q = 0.5, 
we find modifications of the stability properties by about 5%. 

Unlike j3 and 9q, the system dimension a (and the power 
exponent in phonon heat capacity) influences the critical val- 
ues of the scaled activation energy E a jTr considerably since 
it is involved in the Zeldovich number according to Eq. (fTBT l. 
At the same time, the critical Zeldovich number changes 
slightly for different values of a. For example, we obtain the 
scaled activation number E a /Tf = 13.7 and the critical Zel- 



dovich number Z = 6.8 for a = 1, which may be compared 
to E a /Tf — 28.2 and Z = 7.0 found for a = 3. Besides, one 
may also expect different types of heat capacity apart from the 
phonon type, which may modify the power law Eq. (J8j and 
the respective deflagration stability propertie s 9 ' 36,37 . 



VI. SUMMARY 

In this paper we have obtained ID instability of magnetic 
deflagration in a medium of molecular magnets. The main 
parameter of the problem is the Zeldovich number, which 
represents the activation energy of the system (in tempera- 
ture units) scaled by the temperature at the deflagration front 
with a numerical factor depending on the type of heat capac- 
ity. We have demonstrated that the deflagration front becomes 
unstable when the Zeldovich number exceeds a certain criti- 
cal value Z ss 7. We have obtained the analytical scaling for 
the instability parameters at the linear stage within the model 
of an infinitely thin zone of Zeeman energy release. We have 
also solved the eigenvalue stability problem numerically tak- 
ing into account the internal structure of the deflagration front. 
The numerical solution determines the experimental param- 
eters necessary to observe the instability. Besides, we have 
performed direct numerical simulations, which demonstrated 
that the instability leads to a pulsating regime of magnetic de- 
flagration at the nonlinear stage. 
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